#ifndef _QGCM_PSYSTEM_SRC_
#define _QGCM_PSYSTEM_SRC_
#include "psystem.h"
#include "misclib.h"
#include "model.h"
#include "ssystem.h"
//------------------------------------------------------------------------------------------/
psystem::psystem(){
  this->tindex = -1;
  this->mo = NULL;
  this->numtsteps = 0;
  this->model_state.time_initialized = false;
  this->model_state.IC_initialized = false;
};
/* - - - - - - - - - - -- - - - - - - - - - - - - - */
void psystem::init(model_class * model){
  this->mo = model;
  this->numfields = model->get_global_param_by_name(gdb_str_NUMFIELDS).getiValue();
  // set up initial conditions
  this->init_state();
  for(size_t i = 0; i < mean_field_derivatives.size(); i++){
    if( mean_field_derivatives[i].GetPtr() == NULL ){
      throw cgcmexception("Failed to define derivative for field '" + 
          mean_field_names[i] + string("'."), -1 );
    }
  } 
}
/* - - - - - - - - - - -- - - - - - - - - - - - - - */
void psystem::init_state(){
  int iter = 0;
  this->tindex = 0;
  try{
    // initialize the initial fields using parser expressions
    for(iter = 0; iter < static_cast<int>(mean_fields.size()); iter++){
      (mean_fields[iter]) = mean_field_initializers[iter]->Eval();
    }
    this->mean_field_temp.resize( this->mean_fields.size(), 0 );
    this->mean_field_deriv_temp.resize( this->mean_field_derivatives.size(), 0 );
    this->Eval(0, false); // evaluate the first timestep (done again in stepforward, but here
    // it catches an unsuspecting errors in the parsing)
  }catch(const ParserError & e){
    string errormsg(e.what());
    throw cgcmexception(errormsg, -11);
  }
  this->model_state.IC_initialized = true;
}
/* - - - - - - - - - - -- - - - - - - - - - - - - - */
void psystem::addfluxfield(const string & name, INDEX i, INDEX j, const int numfields){
  if(i < 1 || i > numfields || j < 1 || j > numfields){
    fprintf(stderr, "numfields given is %i.\n", numfields);
    throw cgcmexception("Attempted to add s-system flux field with indices that are out of bounds.\n", -1);
    return;
  }
  Value * nval = new Value(static_cast<float_type>(0.0));
  this->sfields.push_back( nval ); // nval is implicitly converted to type Variable
  this->sfield_names.push_back(name);
  rank2 r(i-1, j-1 ); 
  /*r.i = i - static_cast<INDEX>(1); 
    r.j = j - static_cast<INDEX>(1);*/
  this->sfield_indices.push_back(r);
}
/* - - - - - - - - - - -- - - - - - - - - - - - - - */
psystem::~psystem(){
};
/* - - - - - - - - - - -- - - - - - - - - - - - - - */
const unsigned int psystem::getvariablesize() const{
  return this->variables.size();
}
/* - - - - - - - - - - -- - - - - - - - - - - - - - */
const string & psystem::getvariablename(const unsigned int i) const{
  if(i >= 0 && i < variables.size() ){
    return this->variable_names[i];
  }else{
    throw cgcmexception("Fatal error in function 'psystem::getvariablename'. Attempted to access array out of bounds.", -20);
  }
}
/* - - - - - - - - - - -- - - - - - - - - - - - - - */
bool psystem::doesvariableexist(const string & name) const{
  for(unsigned int i = 0; i < this->variables.size();i++){
    if(name == variable_names[i]) return true;
  }
  return false;
}
/* - - - - - - - - - - -- - - - - - - - - - - - - - */
const vector< vector< Value > > & psystem::getvariablehistory() const{
	return this->stored_variables;
}
/* - - - - - - - - - - -- - - - - - - - - - - - - - */
const vector< vector< Value > > & psystem::getmeanfieldhistory() const{
	return this->stored_meanfields;
}

/* - - - - - - - - - - -- - - - - - - - - - - - - - */
const Variable & psystem::getvariablebyindex(const unsigned int i) const{
  if(i >= 0 && i < variables.size() ){
    return variables[i];
  }else{
    fatal_stop("Could not find requested variable.");
    return variables[0];
  }
}
/* - - - - - - - - - - -- - - - - - - - - - - - - - */
int psystem::getfuncbyname(const string & fname) const{
  for(unsigned int i = 0; i < this->function_names.size(); i++){
    if(function_names[i] == fname){ return i; }
  }
  return -1;
}
/* - - - - - - - - - - -- - - - - - - - - - - - - - */
int psystem::getsfieldbyname(const string & sname) const{
  for(unsigned int i = 0; i < this->sfield_names.size(); i++){
    if(this->sfield_names[i] == sname){ return i; }
  }
  return -1;
}
/* - - - - - - - - - - -- - - - - - - - - - - - - - */
int psystem::getfieldbyname(const string & sname) const{
  for(unsigned int i = 0; i < this->mean_field_names.size(); i++){
    if(mean_field_names[i] == sname){ return i; }
  }
  return -1;
}
/* - - - - - - - - - - -- - - - - - - - - - - - - - */
void psystem::init_tindex(const unsigned int & ntsteps){
  if( ! this->model_state.time_initialized ){
    this -> tindex = 0;
  }
  this -> numtsteps = ntsteps;
  assert(ntsteps >= 0);
  Value dNaN( static_cast<float_type>( sqrt(-1) ) );
  vector<Value> ivec;
  if(this->variables.size() > 0){
    ivec.resize(this->variables.size(), dNaN );
    stored_variables.resize(ntsteps+1, ivec); // all other fields are dependent on these and the Ssystem
  }
  // resize all the storage arrays
  if(this->mean_fields.size() > 0){
    ivec.resize(this->mean_fields.size(), dNaN);
    stored_meanfields.resize(ntsteps+1, ivec);
  }
  this->model_state.time_initialized = true;
}
/* - - - - - - - - - - -- - - - - - - - - - - - - - */
void psystem::save_state(int tstep){
  assert( this->model_state.time_initialized );
  rassert( tstep < static_cast<int>(numtsteps) );
  if( variables.size() > 0 ){
    vector<Value> & data = stored_variables[tstep];
    for(size_t i = 0; i < variables.size(); i++){
      data[i] = variables[i].GetFloat();
    }
  }
  if( stored_meanfields.size() > 0 ){
    vector<Value> & data = stored_meanfields[tstep];
    for(size_t i = 0; i < data.size(); i++){
      data[i] = mean_fields[i].GetFloat();
    };
  }
};
/* - - - - - - - - - - -- - - - - - - - - - - - - - */
bool psystem::save_model_state(NcFile * ncfile) const{
  // start by saving all variables in a time array
  for(size_t iter = 0; iter < this->variables.size(); iter++){
    // start by creating the name for the variable field
    string sname("variable_");
    sname += convertint2string(static_cast<int>(iter));
    const string & vname = this->variable_names[iter];
    // copy to an array
    int tsize = tindex + 1;
    double * data = new double[tsize];
    for(size_t ti = 0; ti < static_cast<size_t>(tsize); ++ti) data[ti] = this->stored_variables[ti][iter].GetFloat(); 
    data[tsize-1] = (this->variables[iter].GetFloat());
    NcDim * tdim = ncfile->get_dim("time"); 
    NcVar * var = ncfile->add_var(sname.c_str(), ncDouble, tdim); 
    var->put( data, tsize, 0, 0, 0, 0 );
    var->add_att("name", vname.c_str());
    var->add_att("module", "psystem");
    var->add_att("type", "variable");
    delete [] data;
  }
  // now save all fields
  for(size_t iter = 0; iter < this->mean_fields.size(); iter++){
    // start by creating the name for the variable field
    string sname("field_");
    sname += convertint2string(static_cast<int>(iter));
    const string & vname = this->mean_field_names[iter];
    // copy to an array
    int tsize = tindex + 1;
    double * data = new double[tsize];
    for(size_t ti = 0; ti < static_cast<size_t>(tsize); ++ti) data[ti]=this->stored_meanfields[ti][iter].GetFloat(); 
    data[tsize-1] = (this->mean_fields[iter]).GetFloat();
    NcDim * tdim = ncfile->get_dim("time"); 
    NcVar * var = ncfile->add_var(sname.c_str(), ncDouble, tdim); 
    var->put( data, tsize, 0, 0, 0, 0 );
    var->add_att("name", vname.c_str());
    var->add_att("module", "psystem");
    var->add_att("type", "field");
    delete [] data;
  }
  // now save function information (this must include it's ParserX expression ...)
  // as well as the an index to the mean field it points to if the function is a derivative
  // tricky stuff
  for(size_t iter = 0; iter < functions.size(); iter++){
    // current value is stored as a scalar, if the function is a derivative that information is stored as an attribute
    string sname("function_");
    sname += convertint2string(static_cast<int>(iter));
    const string & vname = this->function_names[iter];
    const string & expr = this->functions[iter]->GetExprAsString();
    const double value = (this->function_values[iter]).GetFloat(); // last value of the function
    NcVar *var = ncfile->add_var(sname.c_str(), ncDouble, NULL); 
    var->put( &(value), 1, 0,0,0,0 );
    var->add_att("name", vname.c_str());
    var->add_att("module", "psystem");
    var->add_att("expression", expr.c_str());
    // see if it is a derivative
    int dindex = -1;
    const Variable & ptr = function_values[iter];
    for( size_t x = 0; x < mean_field_derivatives.size(); x++){
      if( mean_field_derivatives[x] == ptr ){
        dindex = x;
        break;
      }
    }
    var->add_att("type", "function");
    if(dindex >= 0){
      var->add_att("derivative_of", this->mean_field_names[dindex].c_str());
    }
  }
  // finally save the flux fields
  for(size_t iter = 0; iter < sfields.size(); iter++){
    string sname("fluxfield_");
    sname += convertint2string(static_cast<int>(iter));
    const string & vname = this->sfield_names[iter];
    const double lvalue = (this->sfields[iter]).GetFloat();
    const rank2 & r2 = this->sfield_indices[iter];
    NcVar * var = ncfile->add_var(sname.c_str(), ncDouble);
    var->put( &(lvalue),  1, 0, 0, 0, 0 );
    var->add_att("name", vname.c_str());
    var->add_att("module", "psystem");
    var->add_att("type", "fluxfield");
    var->add_att("field_i", static_cast<int>(r2[0]));
    var->add_att("field_j", static_cast<int>(r2[1]));
  } 
  return true;
};
/* - - - - - - - - - - -- - - - - - - - - - - - - - */
void psystem::adjusttime(const unsigned int & ntimes){
  this->numtsteps = ntimes;
  int new_tindex = static_cast<int>( ntimes ) - 1;
  if(new_tindex <= this->tindex){ // do nothing
    this->tindex = new_tindex;
    return;
  }
  this->tindex = new_tindex;
  // else we need to expand the history fields
  int ntsteps = static_cast<int>(ntimes) - 1;
  assert(ntimes >= 0);
  Value dNaN( static_cast<float_type>( sqrt(-1) ) );
  vector<Value> ivec;
  if(this->variables.size() > 0){
    ivec.resize(this->variables.size(), (dNaN) );
    stored_variables.resize(ntsteps+1, ivec); // all other fields are dependent on these and the Ssystem
  }
  // resize all the storage arrays
  if(this->mean_fields.size() > 0){
    ivec.resize(this->mean_fields.size(), dNaN);
    stored_meanfields.resize(ntsteps+1, ivec);
  }
};
// -------------------------------------------------------------------------------------- //
bool restore_vector_history(NcFile * ncfile, const int ntimes,
    const char * label, 
    const vector<string> & variable_names,
    vector<Variable> & variables,
    vector< vector<Value> > & stored_variables){

  vector<bool> check_list( variables.size(), false ); // all variables must be there
  int iter = 0;
  while ( true ){ 
    string nc_var_name(label); nc_var_name+=string("_")+convertint2string(iter);
    NcVar * var = ncfile->get_var(nc_var_name.c_str());
    if( var == NULL ) break;
    char * var_tmp_name = var->get_att("name")->as_string(0);
    const string var_name(var_tmp_name); delete [] var_tmp_name;
    bool found_one = false;
    for(size_t u = 0; u < variable_names.size(); u++){
      if( variable_names[u] == var_name ){ 
        // copy over the data
        double * data = new double[ntimes];
        var->get(data, ntimes, 0,0,0,0);
        for(int ti = 0; ti < ntimes; ti++){
          vector<Value> & vdata = stored_variables[ti];  
          vdata[u] = data[ti];
        }
        (variables[u]) = data[ntimes-1];
        found_one = true;
        check_list[u] = true;
        break;
      }
    }
    if(! found_one ){ // variable existed in saved state but not in current instantiation
      fprintf(stderr, 
          "Errror: Psystem %s '%s' exists in model saved state but not in current model instantiation.\n", label, var_name.c_str());
      return false;
    }
    ++iter;
  }
  // check all variables to make sure they have been defined
  for(size_t u = 0; u < check_list.size(); u++){
    if( !check_list[u] ){
      fprintf(stderr, "Error: Psystem %s '%s' was not found in model saved state.\n", label, variable_names[u].c_str());
      fprintf(stderr, "Cannot load the state of the P-system because %s is missing from file.\n", label);
      return false;
    }
  }

  return true;
};
// -------------------------------------------------------------------------------------- //
bool psystem::restore_model_state(NcFile * ncfile){
  this->init_tindex( 1 );
  NcError error_status( NcError::silent_nonfatal );
  // get the number of times
  NcDim * tdim = ncfile->get_dim("time");
  if(tdim == NULL){
    fprintf(stderr, "Failed to retrieve 'time' dimension. File corrupt?\n"); return false;
  } 
  const int ntimes = static_cast<int>(tdim -> size());
  if(ntimes <= 0){ // bad value
    fprintf(stderr, "Bad 'time' dimension. File corrupt?\n"); return false;
  }
  this->adjusttime( static_cast<unsigned int>( ntimes ) ); // adjust to ntimes numtimesteps
  // now load variables and save their history
  /* - - - - - - - - - - -- - - - - - - - - - - - - - */
  if(!restore_vector_history(ncfile, ntimes, "variable",this->variable_names,this->variables,stored_variables)) return false;
  /* - - - - - - - - - - -- - - - - - - - - - - - - - */
  // load all fields and save their history
  if(!restore_vector_history(ncfile, ntimes, "field", this->mean_field_names, this->mean_fields, stored_meanfields) ) return false;
  /* - - - - - - - - - - -- - - - - - - - - - - - - - */
  // now we need to restore all functions 
  int iter = 0; // all functions in the saved state must exist in current model instantiation
  while(true){
    string nc_var_name("function_"); nc_var_name+=convertint2string(iter);
    NcVar * func = ncfile->get_var(nc_var_name.c_str());
    if( func == NULL ) break;
    NcAtt * afuncname = func->get_att("name");
    char * a_tmp_func_name = afuncname->as_string(0);
    const string funcname(a_tmp_func_name); delete [] a_tmp_func_name;
    // else look up function in current model instantiation
    bool found_one = false;
    for(size_t u = 0; u < function_names.size(); u++){
      if( function_names[u] == funcname ){
        const string & cexpr = functions[u]->GetExprAsString();
        NcAtt * att_old_expr = func->get_att("expression");
        char * old_expr = att_old_expr->as_string(0);
        if( cexpr != old_expr ){
          fprintf(stderr, "Warning: Function '%s' has a different expression in the current model instantiation\n", function_names[u].c_str());
          fprintf(stderr, "         than it has in the model saved state.\n");
          fprintf(stderr, "         Old value: '%s'\n", old_expr);
          fprintf(stderr, "         New value: '%s'\n", cexpr.c_str());
        }
        delete [] old_expr;
        found_one = true;
        break;
      }
    }
    if( ! found_one ){
      fprintf(stderr, "Error: Function '%s' is present in the model saved state but does not exist in\n", funcname.c_str());
      fprintf(stderr, "       in current instantiation of the model.\n");
      fprintf(stderr, "       Original expression for '%s': '%s'\n", funcname.c_str(), func->get_att("expression")->as_string(0));
      return false;
    }
    iter++;
  } 


  /* - - - - - - - - - - -- - - - - - - - - - - - - - */
  // now we need to load all the derivative information (do correct pointer matching)
  // this requires us to reiterate through the list of functions in the ncfile
  iter = 0;
  vector<bool> deriv_check_list( this->mean_fields.size(), false );
  while(true){
    string nc_var_name("function_"); nc_var_name+=convertint2string(iter);
    NcVar * func = ncfile->get_var(nc_var_name.c_str());
    if( func == NULL ) break;
    NcAtt * is_deriv = func->get_att("derivative_of");
    char * afname = NULL;
    if( is_deriv != NULL ){
      afname = is_deriv->as_string(0);
      bool found_one = false;
      for(size_t u = 0; u < this->mean_fields.size(); u++){
        if( mean_field_names[u] == afname ){
          // we must check that derivative matches are correct
          if( this->mean_field_derivatives[u] != this->function_values[iter] ){
            fprintf(stderr, "Function purpose mismatch.\n");
            fprintf(stderr, "Function '%s' is the derivative of field '%s' in saved state but not in\n",
                function_names[iter].c_str(), afname);
            fprintf(stderr, "   current instantiation of the model.\n");
            return false;
          }
          // this->mean_field_derivatives[u] = this->function_values[func_index]; 
          deriv_check_list[u] = true;
          found_one = true;
          break;
        }  
      }
      if(! found_one ){ // uh oh derivative not found in current model instantiation
        fprintf(stderr, "Found derivative function '%s' in model saved state but not in current instantiation of the model.\n",
            afname);
        return false;
      }
      delete [] afname;
    }
    iter++;
  }
  // make sure all mean fields have a derivative loaded
  for(size_t u = 0; u < deriv_check_list.size(); u++){
    if( ! deriv_check_list[u] ){
      fprintf(stderr, "Error: Mean field '%s' does not have derivative information in the saved model state.\n",
          this->mean_field_names[u].c_str());
    }
  }

  /* - - - - - - - - - - -- - - - - - - - - - - - - - */
  // finally we need to make sure that all flux fields from the saved state
  // are present and have the same name as in the current model instantiation
  iter = 0;
  while(true){
    string nc_var_name("fluxfield_"); nc_var_name += convertint2string(iter);
    NcVar * flux_var = ncfile->get_var(nc_var_name.c_str());
    if( flux_var == NULL ) break; // done iterating through saved state flux fields
    NcAtt * a_name = flux_var->get_att("name");
    char * tmp_name = a_name->as_string(0);
    const string fname(tmp_name); delete [] tmp_name;
    bool found_one = false;
    const int a_field_i = flux_var->get_att("field_i")->as_int(0);
    const int a_field_j = flux_var->get_att("field_j")->as_int(0);
    for(size_t u = 0; u < this->sfield_names.size(); u++){
      if( fname == sfield_names[u] ){
        // check field index consistency (must be within range, warning if they don't agree with previous version)
        if( a_field_i < 0 || a_field_i > this->numfields || a_field_j < 0 || a_field_j > this->numfields){
          fprintf(stderr, "Error: Flux '%s' has indices that lie oustide allowed range for the current instantiation of the model.\n",
              fname.c_str());
          fprintf(stderr, "       Please define flux function '%s' with fields i=%i\tj=%i.\n",
              fname.c_str(), a_field_i+1, a_field_j+1);
          return false;
        }
        if( a_field_i != static_cast<int>(sfield_indices[u][0]) || a_field_j != static_cast<int>(sfield_indices[u][1]) ){
          fprintf(stderr, "Warning: Flux '%s' has indices in the current instantiation of the model that\n", fname.c_str());
          fprintf(stderr, "         are different from the model saved state.\n");
          fprintf(stderr, "   Saved state:  i=%i\tj=%i.\n", a_field_i+1, a_field_j+1 );
          fprintf(stderr, "   Curr. state:  i=%i\tj=%i.\n", static_cast<int>(sfield_indices[u][0])+1, 
              static_cast<int>(sfield_indices[u][1])+1);
        }
        found_one = true;
        break;
      }
    }  
    if(!found_one){
      fprintf(stderr, "Error: Flux '%s' was found in model saved state but not found in current instantiation of the model.\n",
          fname.c_str());
      fprintf(stderr, "       Please define flux function '%s' with fields i=%i\tj=%i.\n",
          fname.c_str(), a_field_i+1, a_field_j+1);
      return false;
    }
    iter++;
  }
  this->Eval(this->tindex, false);
  return true;   
};
// -------------------------------------------------------------------------------------- //
int psystem::getnumtsteps() const{
  return static_cast<int>( this->numtsteps );
}
// -------------------------------------------------------------------------------------- //
const int psystem::get_tindex() const{
  return this->tindex;
}
// -------------------------------------------------------------------------------------- //
void psystem::rewindtime(const unsigned int & numsteps){
  if(numsteps > static_cast<unsigned int>(this->tindex)){
    fatal_stop("Cannot rewind past t = 0.");
  }
  unsigned int ctindex = this->tindex;
  this->tindex = ctindex - numsteps;
  // go back and restore the data from that time
  for(unsigned int i = 0; i < mean_fields.size(); i++){
    (mean_fields[i]) = this->stored_meanfields[this->tindex][i];
  }
  for(unsigned int i = 0; i < variables.size(); i++){
    (variables[i]) = this->stored_variables[this->tindex][i];
  }
  ssystem * ssys = &(mo->getSsystem());
  ssys -> rewindtime(numsteps);
};
// -------------------------------------------------------------------------------------- //
void psystem::assignparserdata(ParserX & P) const{
  for(size_t i = 0; i < this->getmeanfieldsize(); i++){
    P.DefineVar( this->getmeanfieldname(i).c_str(), this->getmeanfield(i) );
    // P.DefineVar(this->getmeanfieldname(i).c_str(), const_cast<Variable*>(this->getmeanfield(i)));
    // Note that casting away constness, which is usually frowned upon, is OK here since
    // these variable's values should not be modified within the namelist anyhow. The removal of the
    // constness is necessary to assign the address of the double to a parser variable. 
  }
  // now assign psystem functions as variables
  for(size_t i = 0; i < this->getfunctionssize(); i++){
    P.DefineVar(this->getfunctionname(i).c_str(), this->getfunctionvalue(i));
  }
  // now assign s-system fields as variables 
  for(size_t i = 0; i < this->getsmeanfieldsize(); i++){
    P.DefineVar(this->getsmeanfieldname(i).c_str(),this->getsmeanfield(i));
  }
  // now assign all P-system variables
  for(size_t i = 0; i < this->getvariablesize(); i++){
    P.DefineVar(this->getvariablename(i), this->getvariablebyindex(i)  );
  }
};
// -------------------------------------------------------------------------------------- //
void psystem::Eval(int ti, bool provisional){
  // get s-system fluxes
  ssystem & Ssys = mo->smodule;
  for(size_t xx = 0; xx < this->sfields.size(); xx++){
    const rank2 & r = this->sfield_indices[xx];
    (sfields[xx]) = Ssys.getflux(r[0],r[1], ti, provisional); 
  }
  // evaluate all dependencies (functions, etc);
  for(size_t i = 0; i < functions.size(); i++){ 
    function_values[i]  = functions[i]->Eval(); // (includes derivatives ... )
  }
};

// -------------------------------------------------------------------------------------- //
void psystem::stepforward(REAL timestep, bool nosys, bool nonlinear){
  assert(this->model_state.IC_initialized);
  assert(this->model_state.time_initialized);
  ssystem & Ssys = (mo->getSsystem());
  vector<Value> & mean_fields_0 = this->mean_field_temp;
  vector<Value> & deriv_funcs_0 = this->mean_field_deriv_temp;
  this->save_state( this->tindex );  // store the model state for the current time step
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  // we need to store in order to use them later in the predictor corrector
  // the current value of the mean fields first
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  for(unsigned int i = 0; i < mean_fields.size(); i++){
    mean_fields_0[i] = (mean_fields[i]); // store the initial value for use with Modified Euler later 
    deriv_funcs_0[i] = (mean_field_derivatives[i]); // save all derivatives
  }
  if(!nosys){ //  march the S-system forward to the half (provisional timestep) 
    Ssys.stepforward( static_cast<REAL>(timestep), true); 
  } 
  // now march all P-system fields forward as a provisional step
  for(unsigned int i = 0; i < mean_fields.size(); i++){
    // *(this->mean_fields[i]) = mean_fields_0[i]+*(this->mean_field_derivatives[i])*timestep;   
    REAL old_val = mean_fields_0[i].GetFloat();
    REAL der_val = mean_field_derivatives[i].GetFloat();
    this->mean_fields[i] = (old_val + der_val * timestep) ; 
  } // values stored are now the provisional mean_field values ..
  //  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  // now re-evaluate all dependencies at the intermediate tstep
  // so as the march the full timestep
  //  - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  this->Eval(this->tindex, true); // evalue at half tstep

  if(!nosys){
    Ssys.stepforward(timestep, false); // march the S-system the full timestep
  }
  // store the mean derivatives in the derivative fields
  // apply predictor corrector (modified euler method) 
  for(size_t i = 0; i < mean_fields.size(); i++){
    //*(this->mean_field_derivatives[i]) = 0.5*(deriv_funcs_0[i] + *(this->mean_field_derivatives[i]));
    REAL deriv_func = deriv_funcs_0[i].GetFloat();
    REAL mf_deriv = mean_field_derivatives[i].GetFloat();
    this->mean_field_derivatives[i] = ( 0.5*(deriv_func+ mf_deriv ) );
    //*(this->mean_fields[i]) = mean_fields_0[i] + *(this->mean_field_derivatives[i]) * timestep;
    REAL mf = mean_fields_0[i].GetFloat();
    this->mean_fields[i] = ( mf + mf_deriv * timestep );
  }
  // done stepping to the next time step
  this->tindex ++;
  this->Eval(this->tindex, false); // evalue at the next timestep (not provisional)
}; 
/* - - - - - - - - - - -- - - - - - - - - - - - - - */
const string & psystem::getfunctionexpr(unsigned int index) const{
  if(index >= 0 && index < this->function_names.size()){
    return this->functions[index]->GetExprAsString();
  }else{
    throw cgcmexception("Accessed array out of bounds in psystem::getfunctionexpr.", -1);
    return this->functions[0]->GetExprAsString();
  }
};
/* - - - - - - - - - - -- - - - - - - - - - - - - - */
const string & psystem::getfunctionname(unsigned int index) const{
  if(index >= 0 && index < this->function_names.size()){
    return this->function_names[index];
  }else{
    throw cgcmexception("Accessed array out of bounds in psystem::getfunctionname.", -1);
    return this->function_names[0];
  }
}
/* - - - - - - - - - - -- - - - - - - - - - - - - - */
const unsigned int psystem::getfunctionssize() const{
  return this->functions.size();
}
/* - - - - - - - - - - -- - - - - - - - - - - - - - */
const Variable & psystem::getfunctionvalue(unsigned int index) const{
  if( index < 0 || index >= this->functions.size()){
    fatal_stop("Could not return requested function value.");
  }
  return function_values[index];
}
/* - - - - - - - - - - -- - - - - - - - - - - - - - */
const string & psystem::getmeanfieldinitexpr(unsigned int index) const{
  if(index >= 0 && index < this->mean_field_names.size()){
    return mean_field_initializers[index]->GetExprAsString();
  }else{
    throw cgcmexception("Accessed array out of bounds in psystem::getmeanfieldname.", -1);
    return mean_field_initializers[0]->GetExprAsString();
  }
};
/* - - - - - - - - - - -- - - - - - - - - - - - - - */
const string & psystem::getmeanfieldname(unsigned int index) const{
  if(index >= 0 && index < this->mean_field_names.size()){
    return this->mean_field_names[index];
  }else{
    throw cgcmexception("Accessed array out of bounds in psystem::getmeanfieldname.", -1);
    return this->mean_field_names[0];
  }
}
/* - - - - - - - - - - -- - - - - - - - - - - - - - */
const string & psystem::getsmeanfieldname(unsigned int index) const{
  if(index >= 0 && index < this->sfield_names.size()){
    return this->sfield_names[index];
  }else{
    throw cgcmexception("Accessed array out of bounds in psystem::getsmeanfieldname.", -1);
    return this->sfield_names[0];
  }
}
/* - - - - - - - - - - -- - - - - - - - - - - - - - */
const rank2 & psystem::getsmeanfieldindices(unsigned int index) const{
  if(index >= 0 && index < this->sfield_names.size()){
    return this->sfield_indices[index];
  }else{
    throw cgcmexception("Accessed array out of bounds in psystem::getsmeanfieldname.", -1);
    return this->sfield_indices[0];
  }
}
/* - - - - - - - - - - -- - - - - - - - - - - - - - */
const unsigned int psystem::getmeanfieldsize() const{
  return this->mean_fields.size();
}
/* - - - - - - - - - - -- - - - - - - - - - - - - - */
const unsigned int psystem::getsmeanfieldsize() const{
  return this->sfields.size();
}
/* - - - - - - - - - - -- - - - - - - - - - - - - - */
const Variable & psystem::getsmeanfield(unsigned int index) const{
  if(index >= 0 && index < this->sfields.size()){
    return this->sfields[index];
  }else{
    fatal_stop("Could not get requested smeanfield.");
  }
  return this->sfields[0];
}
/* - - - - - - - - - - -- - - - - - - - - - - - - - */
const Variable & psystem::getmeanfieldderivative(unsigned int index) const{
  if(index >= 0 && index < this->mean_fields.size()){
    return this->mean_field_derivatives[index];
  }
  fatal_stop("Could not get requested mean field derivative.");
  return mean_field_derivatives[0];
};
/* - - - - - - - - - - -- - - - - - - - - - - - - - */
const Variable & psystem::getmeanfield(unsigned int index) const{
  if(index >= 0 && index < this->mean_fields.size()){
    return this->mean_fields[index];
  }
  fatal_stop("Could not get requested mean field.");
  return mean_fields[0];
}
#endif
